---
output: 
  bookdown::pdf_document2:
    fig_caption: yes
    keep_tex: no
    toc: no
    number_sections: yes
    latex_engine: xelatex
    # pandoc_args: --lua-filter=multiple-bibliographies.lua
title: |
  | Revisiting the Evidence on 
  | Thermostatic Response to Democratic Change:
  | Degrees of Democratic Support or 
  | Researcher Degrees of Freedom?
thanks: "Corresponding author: [yhcasstai@psu.edu](mailto:yhcasstai@psu.edu). Current version: `r format(Sys.time(), '%B %d, %Y')`.  Replication materials and complete revision history may be found at [https://github.com/fsolt/dem_mood](https://github.com/fsolt/dem_mood)."
anonymous: true  
keywords: "Democratic support, democracy, measurement, measurement coding"
tables: true # enable longtable and booktabs
citation_package: natbib
fontsize: 12pt
indent: true
linestretch: 1.5 # double spacing using linestretch 1.5
bibliography: dem-mood-text.bib
# csl: https://raw.githubusercontent.com/citation-style-language/styles/master/american-political-science-association.csl
biblio-style: apsr
citecolor: black
linkcolor: black
endnote: no
header-includes:
      - \usepackage{array}
      - \usepackage{caption}
      - \usepackage{graphicx}
      - \usepackage{siunitx}
      - \usepackage{colortbl}
      - \usepackage{multirow}
      - \usepackage{hhline}
      - \usepackage{calc}
      - \usepackage{tabularx}
      - \usepackage{threeparttable}
      - \usepackage{wrapfig}
      - \renewcommand{\topfraction}{.85}
      - \renewcommand{\bottomfraction}{.7}
      - \renewcommand{\textfraction}{.15}
      - \renewcommand{\floatpagefraction}{.66}
      - \usepackage{pdflscape} #\usepackage{lscape} better for printing, page displayed vertically, content in landscape mode, \usepackage{pdflscape} better for screen, page displayed horizontally, content in landscape mode
      - \newcommand{\blandscape}{\begin{landscape}}
      - \newcommand{\elandscape}{\end{landscape}}
      - \setcounter{topnumber}{3}
      - \setcounter{bottomnumber}{3}
      - \setcounter{totalnumber}{4}
      - \usepackage{float}
---



# Authors {-}
-   Yue Hu, ORCID: <https://orcid.org/0000-0002-2829-3971>, Associate Professor, Department of Political Science, Tsinghua University, [yuehu\@tsinghua.edu.cn](mailto:yuehu@tsinghua.edu.cn){.email}
-   Yuehong Cassandra Tai, ORCID: <https://orcid.org/0000-0001-7303-7443>, Postdoctoral Fellow, Center for Social Data Analytics, Pennsylvania State University, [yhcasstai\@psu.edu](mailto:yhcasstai@psu.edu){.email}
-   Frederick Solt, ORCID: <https://orcid.org/0000-0002-3154-6132>, Associate Professor, Department of Political Science, University of Iowa, [frederick-solt\@uiowa.edu](mailto:frederick-solt@uiowa.edu){.email}

\pagebreak
# Abstract {-}

Prominent recent work argues that support for democracy behaves thermostatically---that democratic erosion boosts democratic support while deepening democracy yields public backlash---and further contends that there is no evidence for the classic argument that democracy itself increases democratic support over time.  Here, we document how these conclusions depend on subtle choices in measurement coding that constitute 'researcher degrees of freedom': analyses employing alternative reasonable choices provide little or no support for the original conclusions. The fragility of the statistical results demonstrates that researcher degrees of freedom in measurement must be taken seriously and that the question of the relationship between democratic institutions and democratic support remains unsettled.

\pagebreak

\captionsetup[figure]{list=no}
```{r init, include=FALSE}

options(tinytex.verbose = TRUE)
chooseCRANmirror(graphics = FALSE, ind = 1)

if (!require(pacman))
  install.packages("pacman")
library(pacman)
if (!require(rmarkdown))
  install.packages("rmarkdown")
if (!require(bookdown))
  install.packages("bookdown")
if (!require(knitr))
  install.packages("knitr")

knitr::opts_chunk$set(
  echo = FALSE,
  message = FALSE,
  warning = FALSE,
  cache = TRUE,
  dpi = 300
)

if (!require(cmdstanr)) {
  install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/",
                                         getOption("repos")))
  library(cmdstanr)
  install_cmdstan() # C++ toolchain required; see https://mc-stan.org/cmdstanr/articles/cmdstanr.html
  cmdstanr::check_cmdstan_toolchain(fix = TRUE)
}

if (!require(loo))
  install.packages("loo")
if (!require(BiocManager))
  install.packages("BiocManager")
if (!require(posterior))
  install.packages("posterior")
if (!require(rstan))
  install.packages("rstan")
if (!require(insight))
  install.packages("insight")


p_install(janitor, force = FALSE)
p_install_gh(c("fsolt/DCPOtools"), force = FALSE)

# load all the packages you will use below
p_load(
  # analysis
  cmdstanr,
  plm,
  osfr,
  
  # presentation
  gridExtra,
  modelsummary,
  dotwhisker,
  latex2exp,
  ggh4x,
  
  # data wrangling
  DCPOtools,
  janitor,
  countrycode,
  here,
  broom,
  tidyverse,
  furrr,
  glue)

# Functions preload
set.seed(313)

## data set-up functions
# dichotomize data according to coding rule and missing-data treatment
dichotomize <- function(df, 
                        by = c("highest", "median", "orig", "lowest"),
                        miss = c("mar", "lower", "upper", "theory")) {
  if (miss == "mar") {
    df <- df %>% 
      filter(!r == -1) # missing at random so omitted
  } else if (miss == "upper") { # all missing recoded as supporting democracy 
    df <- df %>% 
      mutate(r = case_when(r == -1 ~ 999,
                           TRUE ~ r))
  } else if (miss == "theory") {
    df <- df %>% 
      mutate(r = case_when(r == -1 & dem == 0 ~ 999, # code non-respondents in 
                           TRUE ~ r))  # autocracies as supporting democracy
  } # else miss == "lower" & all missing answers remain coded as unsupportive
  
  if (by == "highest") {
    df1 <- df %>%
      group_by(country, year, item) %>% 
      summarize(x = sum(ifelse(r >= highest_r, n, 0)),
                samp = sum(n),
                survey = first(survey),
                .groups = "drop_last") %>%
      ungroup() %>%
      filter(x > 0)
  } else if (by == "median") {
    df1 <- df %>% 
      group_by(country, year, item) %>% 
      summarize(x = sum(ifelse(r > median_r, n, 0)),
                samp = sum(n),
                survey = first(survey),
                .groups = "drop_last") %>%
      ungroup() %>%
      filter(x > 0)
  } else if (by == "orig") {
    df1 <- df %>% 
      group_by(country, year, item) %>% 
      summarize(x = sum(ifelse(r > orig_r, n, 0)),
                samp = sum(n),
                survey = first(survey),
                .groups = "drop_last") %>%
      ungroup() %>%
      filter(x > 0)
  } else if (by == "lowest") {
    df1 <- df %>% 
      group_by(country, year, item) %>% 
      summarize(x = sum(ifelse(r > lowest_r, n, 0)),
                samp = sum(n),
                survey = first(survey),
                .groups = "drop_last") %>%
      ungroup() %>%
      filter(x > 0)
  }
  
  df2 <- df1 %>% 
    mutate(by = by,
           miss = miss)
  return(df2)
}

format_claassen <- function(claassen_data) {
  claassen_stan <- map(claassen_data, function(dat) {
    one_stan_input <- list(  N       = nrow(dat),
                             K       = dplyr::n_distinct(dat$item),
                             T       = max(dat$year) - min(dat$year) + 1,
                             J       = dplyr::n_distinct(dat$country),
                             P       = max(as.numeric(as.factor(paste(dat$country, dat$item)))),
                             jj      = as.numeric(as.factor(dat$country)),
                             kk      = as.numeric(as.factor(dat$item)),
                             tt      = dat$year - min(dat$year) + 1,
                             pp      = as.numeric(as.factor(paste(dat$country, dat$item))),
                             x       = round(dat$x),
                             samp    = round(dat$samp),
                             data    = dat)
  })
  return(claassen_stan)
}

## analyzing and presenting functions
# Variance Calculation ------------------------------------------------------------------

## Beck-Katz panel-corrected standard errors
vcovHC_se <-  function(x) {
  plm::vcovHC(x, method="arellano", cluster="group") %>%  #default setting
    diag() %>% 
    sqrt()
}

# Tabulation -----------------------------------------------------------------------
na_types_dict <- list("r" = NA_real_,
                      "i" = rlang::na_int,
                      "c" = NA_character_,
                      "l" = rlang::na_lgl)

# A function that converts a string to a vector of NA types.
# e.g. "rri" -> c(NA_real_, NA_real_, rlang::na_int)
parse_na_types <- function(s) {
  
  positions <- purrr::map(
    stringr::str_split(s, pattern = ""), 
    match,
    table = names(na_types_dict)
  ) %>%
    unlist()
  
  na_types_dict[positions] %>%
    unlist() %>%
    unname()
}

# A function that, given named arguments, will make a one-row
# tibble, switching out NULLs for the appropriate NA type.
as_glance_tibble <- function(..., na_types) {
  
  cols <- list(...)
  
  if (length(cols) != stringr::str_length(na_types)) {
    stop(
      "The number of columns provided does not match the number of ",
      "column types provided."
    )
  }
  
  na_types_long <- parse_na_types(na_types)
  
  entries <- purrr::map2(cols, 
                         na_types_long, 
                         function(.x, .y) {if (length(.x) == 0) .y else .x})
  
  tibble::as_tibble_row(entries)
  
}

tidy.pgmm <- function(x,
                      conf.int = FALSE,
                      conf.level = 0.95,
                      ...) {
  result <- summary(x)$coefficients %>%
    tibble::as_tibble(rownames = "term") %>%
    dplyr::rename(
      estimate = Estimate,
      std.error = `Std. Error`,
      statistic = `z-value`,
      p.value = `Pr(>|z|)`
    )
  
  if (conf.int) {
    ci <- confint(x, level = conf.level) %>% 
      as.data.frame() %>% 
      rownames_to_column(var = "term") %>% 
      dplyr::rename(
        conf.low = `2.5 %`,
        conf.high = `97.5 %`
      )
    result <- dplyr::left_join(result, ci, by = "term")
  }
  
  result
}

glance.plm <- function(x, ...) {
  s <- summary(x)
  as_glance_tibble(
    nobs = stats::nobs(x),
    n.country = pdim(x)$nT$n,
    na_types = "ii"
  )
}

glance.pgmm <- function(x, ...) {
  s <- summary(x)
  as_glance_tibble(nobs = stats::nobs(x),
                   n.country = pdim(x)$nT$n,
                   n.inst = dim(x$W[[1]])[2],
                   na_types = "iii"
  )
}

theme_set(theme_bw())

## dotwhisker::small_multiple() hacks
body(small_multiple)[[19]] <- substitute(
  p <-
    ggplot(
      df,
      aes(
        y = estimate,
        ymin = conf.low,
        ymax = conf.high,
        x = as.factor(model),
        colour = submodel
      )
    ) +
    do.call(geom_pointrange, point_args) +
    ylab("") + xlab("") +
    facet_grid(
      term ~ .,
      scales = "free",
      labeller = label_parsed,
      # enable LaTeX facet labels
      switch = "y"
    ) +             # put facet labels on left
    scale_y_continuous(position = "right") # put axis label on right
)
```

```{r claassen_input, include=FALSE}

# publication data
sd <- read_csv(here::here("data", "supdem raw survey marginals.csv"), col_types = "cdcddcdc") %>% 
  mutate(item_fam = str_extract(tolower(Item), "^[a-z]+"),
         country = countrycode::countrycode(Country, "country.name", "country.name"),
         dataset = "supdem") %>% 
  rename(year = Year, project = Project) %>% 
  DCPOtools::with_min_yrs(2) # Selecting data w. at least two years

# 
if (!file.exists(here::here("data", "claassen_input_raw.csv"))) {
  claassen_input_raw <- DCPOtools:::claassen_setup(
    vars = read_csv(here::here("data-raw", "mood_dem.csv"),
                    col_types = "cccccc")) %>% 
    pluck("lower")
  
  write_csv(claassen_input_raw, 
            file = here::here("data",
                              "claassen_input_raw.csv"))
}

claassen_input_raw0 <- read_csv(here::here("data",
                                           "claassen_input_raw.csv"),
                                col_types = "cdcddcd")

claassen_input_raw1 <- claassen_input_raw0 %>%
  filter(!((
    str_detect(item, "army_wvs") &
      # WVS obs identified as problematic by Claassen
      ((country == "Albania" & year == 1998) |
         (country == "Indonesia" &
            (year == 2001 | year == 2006)) |
         (country == "Iran" & year == 2000) |
         (country == "Pakistan" &
            (year == 1997 | year == 2001)) | # 1996 in Claassen
         (country == "Vietnam" & year == 2001)
      )
  ) |
    (
      str_detect(item, "strong_wvs") &
        ((country == "Egypt" & year == 2012) |
           (country == "Iran" &
              (year == 2000 | year == 2007)) | # 2005 in Claassen
           (country == "India") |
           (country == "Pakistan" &
              (year == 1997 | year == 2001)) | # 1996 in Claassen
           (country == "Kyrgyzstan" &
              (year == 2003 | year == 2011)) |
           (country == "Romania" &
              (year == 1998 | year == 2005 | year == 2012)) |
           (country == "Vietnam" & year == 2001)
        )
    ) |
    (
      country %in% c(
        "Puerto Rico",
        "Northern Ireland",
        "SrpSka Republic",
        "Hong Kong SAR China"
      )
    ))) %>%
  with_min_yrs(2)

claassen_input0 <- DCPOtools::format_claassen(claassen_input_raw1)

cri <- claassen_input0$data %>% 
  mutate(p_dcpo = str_extract(survey, "^[a-z]+"), 
         project = case_when(p_dcpo == "afrob" ~ "afb",
                             p_dcpo == "amb" ~ "lapop",
                             p_dcpo == "arabb" ~ "arb",
                             p_dcpo == "asiab" ~ "asb",
                             p_dcpo == "asianb" ~ "asnb",
                             p_dcpo == "neb" ~ "ndb",
                             p_dcpo == "sasianb" ~ "sab",
                             TRUE ~ p_dcpo),
         item_fam = str_extract(item, "^[a-z]+"),
         item_fam = if_else(item_fam == "election", "elec", item_fam),
         dataset = "cri")

supdem_cri <- full_join(sd, cri, by = c("country", "year", "item_fam", "project"))

# While DCPO assigns actual year of survey, Claassen (2020) often uses nominal year
# of wave, so some corrections are required to match observations (~8% of obs)
no_problems <- inner_join(sd %>% select(-dataset), 
                          cri %>% select(-dataset),
                          by = c("country", "year", "item_fam", "project")) # 3400 obs

needed <- anti_join(sd %>% select(-dataset),
                    cri %>% select(-dataset))                   # 316 obs

available <- anti_join(cri %>% select(-dataset),
                       sd %>% select(-dataset))                # 1336 obs

year_fixes <- left_join(needed,    # 304 obs
                        available,
                        by = c("country", "project", "item_fam")) %>% 
  mutate(diff = year.x - year.y) %>% 
  group_by(country, project, item_fam, year.x) %>% 
  mutate(closest_to_claassen = min(abs(diff))) %>% 
  ungroup() %>% 
  group_by(country, project, item_fam, year.y) %>% 
  mutate(closest_to_dcpo = min(abs(diff))) %>% 
  ungroup() %>% 
  filter(closest_to_claassen == abs(diff) & closest_to_dcpo == abs(diff) & abs(diff) <= 3) %>% 
  filter(!(country == "Egypt" & year.x == 2014 & survey == "afrob5")) %>%  # double matches (it's really afrob6)
  mutate(year = year.y)

claassen_sample <- bind_rows(no_problems, year_fixes) %>% 
  select(country, year, item)

df_apsr <- read_csv(here::here("data", "dem_mood_apsr.csv"),
                    show_col_types = FALSE)

base_input_raw <- claassen_input_raw0 %>%
  right_join(claassen_sample,
             by = join_by(country, year, item)) %>% 
  left_join(df_apsr %>%
              transmute(country = countrycode(Country,
                                              "country.name",
                                              "country.name"),
                        year = Year,
                        dem = Regime_di),
            by = join_by(country, year)) %>% 
  group_by(country, year, item) %>% 
  mutate(highest_r = max(r),
         lowest_r = 1,
         median_r = median(1:first(highest_r)),
         orig_r = case_when((str_detect(survey, "asiab") &
                               item == "evdemoc_asiab") ~ 1,
                            (str_detect(survey, "pew") &
                               item == "impdemoc_pew") ~ 3,
                            TRUE ~ median_r)) %>% 
  ungroup()

by <- c("orig", "highest", "median", "lowest") %>% # coding rules
  set_names()

miss <- c("lower", "upper", "theory", "mar") %>% # missing-data treatments
  set_names()

claassen_input_raw <- map(by,
                          function(b) {
                            map(miss,
                                function(m) dichotomize(base_input_raw, b, m))
                          })

claassen_input <- map(by, function(b) {
  rule <- claassen_input_raw %>% 
    pluck(b) 
  format_claassen(rule)
})

rio::export(claassen_input, 
            file = here("data",
                        "claassen_input.rds"))
```

With democracy under threat in many countries around the world, how the public reacts to democratic change is a crucial question.
A long-established and indeed still-vibrant literature holds that experiencing more democracy boosts democratic support in the public---in short, that democracy creates its own demand [see, e.g., @Lipset1959; @Welzel2013; @Wuttke2022].
In contrast to this classic theory, one prominent recent article, @Claassen2020b, argues that democratic support behaves thermostatically: that increases in democracy yield an authoritarian backlash in the public, while democratic backsliding prompts the public to rally to democracy's cause.

In support of this thermostatic argument, @Claassen2020b offers evidence based on recent advances in modeling public opinion as a latent variable.
Drawing on aggregated responses to more than sixty distinct questions on attitudes toward democracy and its alternatives that were asked in thousands of national surveys, the article estimates a latent variable of democratic support in more than a hundred countries over as long as thirty years.
These data constitute a broader evidentiary base than that employed in any earlier work on this topic.
Analyses of these data support the article's conclusions that democratic support behaves thermostatically in the public, that "increases in democracy dampen public mood, while decreases cheer it," and that democracy "does not appear to create its own demand" as the classic theory has it [@Claassen2020b, 51].

But, inevitably, the aggregated survey data, the resulting estimates of democratic support, and ultimately the results of the article's analyses embody a particular set of choices.
These choices, known as "researcher degrees of freedom," have attracted growing attention among political scientists in recent years [see, e.g., @Wuttke2019; @Breznau2022].
The concern extends beyond cases of '$p$-hacking,' in which researchers sift through many different options to find a set of choices that yield statistically significant results.
It turns out that no intentional fraud or misconduct is required to render results unreliable.
Researchers may make only a single set of entirely reasonable choices, that is, they may take only a single walk through the metaphorical "garden of forking paths"; find a result that supports their expectations; and "get excited and believe it" without even considering that a different set of entirely reasonable choices would provide different results [@Gelman2014, 464].

Here, we illustrate the importance of taking seriously researcher degrees of freedom by replicating the primary analyses of @Claassen2020b while varying two aspects of the measurement of democratic support.
The first is the manner in which ordinal survey responses are coded as supportive or not supportive of democracy.
The second is the attitude ascribed to survey respondents who did not answer items regarding their support for democracy.
We document a range of reasonable alternatives to the choices made in @Claassen2020b on each of these issues, and demonstrate that analyses employing these alternatives provide little to no support for the conclusions drawn in that article.
The fragility of the published results to alternate reasonable choices demonstrates the importance of taking researcher degree of freedom seriously and that the question of the relationship between democratic institutions and democratic support remains unsettled.


## Researcher Degrees of Freedom and Democratic Support {.unlisted .unnumbered}

The dependent variable considered in @Claassen2020b is the public's support for democracy, a latent variable.
This latent variable was estimated from aggregated responses to many survey questions asked in many countries in many years.
Here we examine the consequences of two choices made early in the process of collecting and aggregating all of those survey data: one regarding how responses were coded as supporting democracy and a second involving how missing responses were treated.

First we consider the question of how to code survey responses.
The latent variable model used in @Claassen2020b requires responses to be dichotomous: either, in this case, democracy supporting or not.
This is straightforward for the few survey items that gave respondents only two options, for example, in the Pew Global Attitudes surveys, "a democratic form of government" or "a leader with a strong hand."
But the vast majority of the survey items employed are ordinal and so require a cutpoint to be chosen to split respondents between those who support democracy and those who do not.

There are at least three reasonable options for recoding these ordinal responses into a dichotomous variable of democratic support.
The first coding rule is the most demanding.
It stipulates that only those who supply only the most democracy-supporting response of those available can be considered unflinching supporters of democracy; anything but the highest response indicates a lack of support.
A second coding rule takes the opposite tack, with all but the lowest response considered as indicating at least some support for democracy and only the lowest value considered as not supportive.
The last splits the difference, considering answers above the median, e.g., the two highest values on a five-point scale, to indicate that the respondent supports democracy, with the median and below being unsupportive.

In the event, @Claassen2020b resorts to a mix of the above three rules.
Most survey items were dichotomized following the 'above the median' rule [see @Claassen2020b, Appendix 1.3]. 
Although not documented in the article's appendix on "Microlevel Coding of Survey Responses," the other two rules were employed as well.
The most demanding rule was applied to the four-point item asked in Pew Global Attitudes surveys, "How important is it to you to live in a country where honest elections are held regularly with a choice of at least two political parties? Is it very important, somewhat important, not too important or not important at all?"
For this question, rather than including respondents who gave both responses above the median---"very important" and "somewhat important"---only those respondents who answered "very important" were entered as supporting democracy.
And for the three-point item in the Asia Barometer asking whether respondents thought "a democratic political system" would be very good, fairly good, or bad for their country, the 'all but the lowest' rule was applied.^[\doublespacing
One charitable reading, offered by @Hu2022, is that these deviations from the 'above the median' rule documented in the paper's appendix are simply unintentional mistakes.
A corrigendum to the article written in response to @Hu2022, however, rejects this interpretation [see @Claassen2023].]
Adding this mixed rule to those described above gives us a total of four coding rules.

The second choice involves how to treat missing responses, that is, when survey respondents either indicate that they "don't know" or simply refuse to answer the question.
Here, too, there is a range of reasonable choices.
We identify four possibilities.
One might reasonably assume that a missing response is equivalent to answering with a lack of support for democracy since the respondent did not provide a democracy-supporting response.
@Claassen2020b [, Appendix 1.3] takes this first approach.
But one might conversely and equally reasonably consider a missing response to suggest _support_ for democracy; the respondent, after all, did not provide a democracy-opposing response.
A third possibility is that non-responses indicate opposition to the current regime. That is, their meaning depends on the context.
In democracies missing responses indicate an unwillingness to supply an honest but socially unacceptable rejection of democracy, but in autocracies missing responses evince support for democracy coupled with a fear of reprisal.
One final possibility is that non-responses occur at random.
In that case, a non-response tells nothing about the respondent's views toward democracy and can simply be dropped from the sample.
There are, then, at least four reasonable treatments of survey non-responses.
Together, the four coding rules and four missing-data treatments described above give us sixteen combinations of plausible researcher choices.


## Consequences for Inference {.unlisted .unnumbered}

```{r cm5_output, eval=FALSE, cache=FALSE, include=FALSE, results=FALSE}
# omit eval=FALSE to re-generate the corrected democracy support series
# not evaluated by default because it takes a long while on most machines
cm5 <- cmdstan_model(here("R", "supdem.stan.mod5.stan"))

warmup <- 1000; sampling <- 200

plan(multisession, workers = 12)
future_options <- furrr_options(seed = 324)

orig_output <- future_map(claassen_input %>% 
                            pluck("orig"),
                          \(ci_df) {
                            cm5$sample(
                              data = ci_df, 
                              max_treedepth = 14,
                              adapt_delta = 0.99,
                              step_size = 0.005,
                              seed = 324, 
                              chains = 3, 
                              parallel_chains = 3,
                              iter_warmup = warmup,
                              iter_sampling = sampling,
                              refresh = warmup/25
                            )
                          },
                          .options = future_options
)

results_path <- here::here(file.path("data", 
                                     "orig"))
dir.create(results_path, 
           showWarnings = FALSE, 
           recursive = TRUE)
map(miss, function(m) {
  dir.create(file.path(results_path, m),
             showWarnings = FALSE, 
             recursive = TRUE)
  mt <- orig_output %>% 
    pluck(m)
  mt$save_data_file(dir = file.path(results_path, m),
                    random = FALSE)
  mt$save_output_files(dir = file.path(results_path, m),
                       random = FALSE)
})

highest_output <- future_map(claassen_input %>% 
                               pluck("highest"),
                             \(ci_df) {
                               cm5$sample(
                                 data = ci_df, 
                                 max_treedepth = 14,
                                 adapt_delta = 0.99,
                                 step_size = 0.005,
                                 seed = 324, 
                                 chains = 3, 
                                 parallel_chains = 3,
                                 iter_warmup = warmup,
                                 iter_sampling = sampling,
                                 refresh = warmup/25
                               )
                             },
                             .options = future_options
)


results_path <- here::here(file.path("data", 
                                     "highest"))
dir.create(results_path, 
           showWarnings = FALSE, 
           recursive = TRUE)
map(miss, function(m) {
  dir.create(file.path(results_path, m),
             showWarnings = FALSE, 
             recursive = TRUE)
  mt <- highest_output %>% 
    pluck(m)
  mt$save_data_file(dir = file.path(results_path, m),
                    random = FALSE)
  mt$save_output_files(dir = file.path(results_path, m),
                       random = FALSE)
})

median_output <- future_map(claassen_input %>% 
                              pluck("median"),
                            \(ci_df) {
                              cm5$sample(
                                data = ci_df, 
                                max_treedepth = 14,
                                adapt_delta = 0.99,
                                step_size = 0.005,
                                seed = 324, 
                                chains = 3, 
                                parallel_chains = 3,
                                iter_warmup = warmup,
                                iter_sampling = sampling,
                                refresh = warmup/25
                              )
                            },
                            .options = future_options
)


results_path <- here::here(file.path("data", 
                                     "median"))
dir.create(results_path, 
           showWarnings = FALSE, 
           recursive = TRUE)
map(miss, function(m) {
  dir.create(file.path(results_path, m),
             showWarnings = FALSE, 
             recursive = TRUE)
  mt <- median_output %>% 
    pluck(m)
  mt$save_data_file(dir = file.path(results_path, m),
                    random = FALSE)
  mt$save_output_files(dir = file.path(results_path, m),
                       random = FALSE)
})

lowest_output <- future_map(claassen_input %>% 
                              pluck("lowest"),
                            \(ci_df) {
                              cm5$sample(
                                data = ci_df, 
                                max_treedepth = 14,
                                adapt_delta = 0.99,
                                step_size = 0.005,
                                seed = 324, 
                                chains = 3, 
                                parallel_chains = 3,
                                iter_warmup = warmup,
                                iter_sampling = sampling,
                                refresh = warmup/25
                              )
                            },
                            .options = future_options
)


results_path <- here::here(file.path("data", 
                                     "lowest"))
dir.create(results_path, 
           showWarnings = FALSE, 
           recursive = TRUE)
map(miss, function(m) {
  dir.create(file.path(results_path, m),
             showWarnings = FALSE, 
             recursive = TRUE)
  mt <- lowest_output %>% 
    pluck(m)
  mt$save_data_file(dir = file.path(results_path, m),
                    random = FALSE)
  mt$save_output_files(dir = file.path(results_path, m),
                       random = FALSE)
})
```

```{r df_apsr_list, eval=FALSE}

for (i in by) {
  for (j in miss) {
    # Construct the path to check
    path_to_check <- file.path(here::here("data", i, j))
    
    if (!dir.exists(path_to_check)) {
      dir.create(path_to_check, recursive = TRUE, showWarnings = TRUE)
      message("Created directory: ", path_to_check)
    }
  
    csv_files <- list.files(path = path_to_check, pattern = "\\.csv$", full.names = TRUE)
    if(length(csv_files) == 0) {
      # If the file does not exist, download from OSF
      message("Downloading: ", path_to_check)
      # Retrieve the OSF node and list files
      osf_retrieve_node("vmxkn") %>%
        osf_ls_files(path = file.path(i, j)) %>%
        osf_download(path = path_to_check)
    } else {
      message("File exists: ", path_to_check)
    }
  }
}

load_results <- function(by) {
 # b_results_path <- file.path("data", by)
  map(miss,
      function(m) {
        bm_results_path <- file.path("data", by, m)
        as_cmdstan_fit(here(bm_results_path,
                            list.files(here(bm_results_path), pattern = "csv$")))
      })
}  

orig_output <- load_results("orig")
highest_output <- load_results("highest")
median_output <- load_results("median")
lowest_output <- load_results("lowest")

summarize_cm5_results <- function(cm5_input,
                                  cm5_output) {
  
miss <- c(lower = "lower", upper = "upper", mar = "mar", theory = "theory")
  
res <- map(miss, function(m) {  
    dat <- cm5_input %>% 
      pluck(m) %>% 
      pluck("data")
    
    kcodes <- dat %>%
      dplyr::transmute(country = country,
                       kk = as.numeric(as.factor(country))) %>% 
      unique()
    
    tcodes <- tibble(year = min(dat$year):max(dat$year),
                     tt = year - min(year) + 1)
    
    ktcodes <- dat %>%
      dplyr::group_by(country) %>%
      dplyr::summarize(first_yr = min(year),
                       last_yr = max(year))
    
    fit <- cm5_output %>% 
      pluck(m)
    
    summary_measures <- c("mean", "median", "sd", "mad", "rhat", "ess_bulk", "ess_tail")
    
    suppressWarnings({
      fit$summary("theta",
                  ~posterior::quantile2(., probs = c(.1, .9)),
                  summary_measures) %>%
        dplyr::mutate(tt = as.numeric(gsub("theta\\[(\\d+),\\d+\\]",
                                           "\\1",
                                           variable)),
                      kk = as.numeric(gsub("theta\\[\\d+,(\\d+)\\]",
                                           "\\1",
                                           variable))) %>%
        dplyr::left_join(kcodes, by = "kk") %>%
        dplyr::left_join(tcodes, by = "tt") %>%
        dplyr::left_join(ktcodes, by = "country") %>%
        dplyr::filter(year >= first_yr) %>% # & year <= last_yr
        dplyr::arrange(kk, tt) %>%
        dplyr::select(country, year,
                      starts_with("me"),
                      sd, mad, starts_with("q"),
                      rhat, starts_with("ess"),
                      variable, kk, tt)
    })
  })
  
  return(res)
}

claassen_input <- rio::import(here("data",
                                   "claassen_input.rds"))

theta_summary <- map(by, function(b) {
  summarize_cm5_results(claassen_input[[b]],
                        get(paste0(b, "_output")))
})


if (!file.exists(here::here("data", "dem_mood_apsr.csv"))) {
  tempfile <- dataverse::get_file_by_doi(filedoi = "doi:10.7910/DVN/FECIO3/ACMGQG") # APSR replication
  
  writeBin(tempfile, here::here("data", "dem_mood_apsr.csv"))
  rm(tempfile)
}

df_apsr <- read_csv(here::here("data", "dem_mood_apsr.csv"),
                    show_col_types = FALSE) %>% 
  mutate(Country = countrycode(Country,
                               "country.name",
                               "country.name"))

df_apsr_list <- map(by,
                    function(b) {
                      map(miss,
                          function(m) {
                            df_apsr %>%
                              left_join(theta_summary[[b]][[m]], 
                                        by = c("Year" = "year",
                                               "Country" = "country")) %>% 
                              mutate(SupDem_trim = mean) %>% 
                              select(Country,
                                     Year,
                                     SupDem_trim,
                                     Libdem_z,
                                     lnGDP_imp) %>% 
                              filter(., complete.cases(.))
                          })
                    })

df_apsr_list[[1]][[1]] <- df_apsr %>% 
  select(Country, Year, SupDem_trim, Libdem_z, lnGDP_imp) %>% 
  filter(., complete.cases(.))

rio::export(df_apsr_list, file = here("data",
                                      "df_apsr_list.rds"))
```

How does the selection among these sixteen garden paths influence our conclusions regarding the public's response to changes in democracy?
To find out, we collected all of the survey data used in @Claassen2020b from their original sources.
We then applied each combination of coding rule and non-response treatment and generated sixteen sets of cross-national time-series estimates of democratic support using the latent-variable model employed in the article.
Finally, we used each of these sets of estimates to replicate the article's analyses.

```{r regression_results}
df_apsr_list <- rio::import(here("data", "df_apsr_list.rds"))

by_names <- paste(c("Mixed",
                    "Above Median",
                    "Only Highest",
                    "All But Lowest"),
                  "Coding Rule")

miss_names <- c("Unsupportive",
                    "Supportive",
                    "Oppositional",
                    "At Random")

m1_1_list <- map(by,
                 function(b) {
                   bm <- map(miss,
                       function(m) {
                         plm(diff(SupDem_trim, lag = 1) ~ 
                               plm::lag(SupDem_trim, 1:2) +
                               diff(Libdem_z, lag = 1) +
                               plm::lag(Libdem_z, 1) +
                               diff(lnGDP_imp, lag = 1) +
                               plm::lag(lnGDP_imp, 1),
                             model = 'pooling',
                             data = pdata.frame(df_apsr_list[[b]][[m]],
                                                index = c("Country", "Year")))
                       })
                   names(bm) <- miss_names
                   return(bm)
                 })

names(m1_1_list) <- by_names

m1_1_tidy <- map(by_names,
                 function(b) {
                   map(miss_names,
                       function(m) {
                         tidy_result <- tidy(m1_1_list[[b]][[m]],
                                             conf.int = TRUE) %>% 
                           mutate(std.error = vcovHC_se(m1_1_list[[b]][[m]]),
                                  model = paste(b, m),
                                  coding = b,
                                  treatment = m)
                       }) %>% 
                     bind_rows()
                 }) %>% 
  bind_rows()
```


```{r regressionPlot, fig.cap= "The Effects of Democracy on the Change in Public Support", fig.width=7, fig.height=9}
txt_caption <- strwrap("Notes: Replications of Claassen (2020, 47), Table 1, Model 1.1.  The mixed coding rule employed in Claassen (2020) along with that work's assumption that non-responses indicate a lack of support for democracy yield a larger negative point estimate of the coefficient for change in liberal democracy than most other combinations and a larger point estimate of the coefficient for the lagged level of liberal democracy than all other combinations.  In error-correction models like these, both coefficients must be interpreted together; see Figure 2.", width = 110) %>% 
  paste0(sep="", collapse="\n")

m1_1_tidy %>%
  filter(str_detect(term, "Libdem")) %>% 
  mutate(term = if_else(term == "diff(Libdem_z, lag = 1)",
                        "Delta~Liberal~Democracy",
                        "Liberal~Democracy[t-1]"),
         model = factor(treatment,
                        levels = miss_names,
                        labels = c("Unsupportive",
                                   "Supportive",
                                   "Oppositional",
                                   "At Random")), 
         submodel = factor(coding,
                           levels = by_names,
                           labels = c("Mixed",
                                      "Above Median",
                                      "Only Highest",
                                      "All But Lowest"))) %>%
  small_multiple(model_order = c("Unsupportive",
                                 "Supportive",
                                 "Oppositional",
                                 "At Random")) +
  geom_hline(yintercept = 0, colour = "grey50", linetype = 2) +
  scale_color_viridis_d(name = "Coding Rule",
                        option = "magma",
                        end = .8) +
  theme(strip.text.y.left = element_text(angle = 0),
        strip.background = element_rect(colour="white", fill="white"),
        legend.position = c(-.01, 1),
        legend.justification = c(1, 1), 
        legend.title.align = 0.5,
        legend.text = element_text(size = 8),
        legend.background = element_rect(color="gray90"),
        legend.spacing = unit(-5, "pt"),
        legend.key.width = unit(4.5, "mm"),
        legend.key.height = unit(7, "mm"),
        plot.title.position = "plot",
        plot.caption.position = "plot",
        plot.caption = element_text(size = 10, hjust = 0)) +
  guides(x = guide_axis(title = "Non-Response Treatment")) +
  ggtitle("Predicting the Change in Public Democratic Support") +
  labs(caption = txt_caption)
```

Figure \@ref(fig:regressionPlot) is a "small multiple" plot [see @SoltHu2015] showing the results of replicating Model 1.1, the principal model of @Claassen2020b [, 47], with each of the sixteen combinations of coding rule and treatment of survey non-response.
In the top panel, the dots represent point estimates for the coefficients for change in liberal democracy; according to the thermostatic theory, these coefficients should be negative.
In the bottom panel, the dots depict point estimates for the coefficients for lagged level of liberal democracy; according to the classic theory, these coefficients should be positive.
In both panels, the whiskers show the associated 95% confidence intervals.
Each coding rule is represented by a different color, while the four non-response treatments are shown in separate grouped estimates from left to right.^[
The full tabular results can be found in the online Supplementary Materials.]

The darkest, left-most coefficients in the left-most group of estimates represent the combination of mixed coding rule and treatment of non-responses as unsupportive of democracy, that is, the combination employed in @Claassen2020b.
It replicates the results reported in that article exactly.
More importantly for our purposes, this combination yields a larger negative point estimate of the coefficient for change in liberal democracy than most other combinations.
It also yields a larger point estimate of the coefficient for the lagged level of liberal democracy than all of the other combinations.
But like the coefficients of constitutive terms of multiplicative interactions [see, e.g., @Brambor2006], the coefficients of error-correction models like Model 1.1 are properly interpreted together rather than separately [see @Williams2012].

```{r simulations}
sim_dem_mood <- function(m1_1_result) {
  mod <- m1_1_result
  n.coef =  length(coef(mod)) # 7
  n.yrs = 231
  n.burn = 200
  pred.data = data.frame(Intercept=1,
                         SupDem_m1=0,
                         SupDem_m2=0, 
                         ChgDem=0,
                         Libdem_m1=c(rep(-.5, n.burn),
                                     rep(.5, n.yrs-n.burn)), 
                         ChglogGDP=0,
                         lnGDP_m1=mean(df_apsr_list %>% 
                                         pluck("median") %>% 
                                         pluck("lower") %>% 
                                         pull(lnGDP_imp), na.rm=TRUE),
                         ChgSup=NA,
                         SupDem=NA)
  
  for(i in 1:(n.yrs-1)) {
    pred.data[i, "ChgDem"] = pred.data[i+1, "Libdem_m1"] - pred.data[i, "Libdem_m1"]
    pred.data[i, "ChgSup"] = sum(pred.data[i, 1:n.coef] * coef(mod))
    pred.data[i, "SupDem"] = pred.data[i, "SupDem_m1"] + pred.data[i, "ChgSup"]
    pred.data[i+1, "SupDem_m1"] = pred.data[i, "SupDem"]
    pred.data[i+1, "SupDem_m2"] = pred.data[i, "SupDem_m1"]
  }
  
  n.sims = 5e4
  mod.sims = MASS::mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovHC(mod))
  yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
  sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
  for(i in 1:(n.yrs-1)) {
    yhat[, i] = mod.sims %*% t(pred.data[i, 1:n.coef])
    sup.hat[, i] = pred.data[i, "SupDem_m1"] + 
      yhat[, i] +
      rnorm(n.sims, 0, sigma(mod))
  }
  
  sup_hat <- sup.hat - apply(sup.hat, 2, mean, na.rm=TRUE)[n.burn-1]
  
  sim_data11 <- tibble(fd = apply(sup_hat, 2, mean, na.rm=TRUE),
                       u95 = apply(sup_hat, 2, quantile, 0.975, na.rm=TRUE),
                       l95 = apply(sup_hat, 2, quantile, 0.025, na.rm=TRUE),
                       model = "Model 1.1",
                       dem = "Liberal Democracy") %>% 
    mutate(year = row_number() - n.burn) %>% 
    filter(year > -5 & !is.na(fd))
  
  return(sim_data11)
}

sim_m1_1_list <- map(by_names,
                     function(b) {
                       map(miss_names,
                           function(m) {
                             sim_dem_mood(m1_1_list[[b]][[m]]) %>% 
                               mutate(by = factor(b,
                                                  levels = by_names,
                                                  labels = c("Mixed",
                                                             "Above Median",
                                                             "Only Highest",
                                                             "All But Lowest")),
                                      miss = factor(m,
                                                    levels = miss_names,
                                                    labels = c("Unsupportive",
                                                               "Supportive",
                                                               "Oppositional",
                                                               "At Random")))
                           }) %>% 
                         bind_rows()
                     }) %>% 
  bind_rows()
```

```{r simulationPlot, fig.cap = "Simulated Effects of Democracy on Changes in Public Democratic Support", fig.width=7, fig.height=9}
txt_caption <- strwrap("Notes: Simulated effects are estimated using coefficients from the models presented in Figure 1. The solid lines indicate the mean simulated effect; the shaded regions indicate the 95% confidence intervals of these effects.", width = 115) %>% 
  paste0(sep="", collapse="\n")

ggplot(data = sim_m1_1_list) +
  geom_hline(aes(yintercept = 0), colour="gray50", linetype="dashed") +
  geom_line(aes(x = year, y = fd)) +
  geom_ribbon(aes(x = year, ymin = l95, ymax = u95, alpha = .6)) +
  facet_grid(rows = vars(miss),
             cols = vars(by),
             switch = "y") +
  scale_y_continuous(position = "right") +
  theme(strip.background = element_rect(colour="white", fill="white"),
        legend.position = "none",
        plot.caption.position = "plot",
        plot.caption = element_text(size = 10, hjust = 0)) +
  ggtitle("Simulating the Effect of a Change in Democracy on Public Democratic Support") +
  guides(x = guide_axis(title = "Year"),
         x.sec = guide_none(title = "Coding Rule"),
         y = guide_axis(title = element_blank()),
         y.sec = guide_none(title = "Non-Response Treatment")) +
  labs(caption = txt_caption)
```

Figure \@ref(fig:simulationPlot) is similar to Claassen's [-@Claassen2020b, 48] Figure 5.
It depicts simulated effects, in differences rather than levels for ease of interpretation, of a one standard deviation increase in democracy on the public's support for democracy using the sixteen sets of regression coefficients presented in Figure \@ref(fig:regressionPlot) with the four different coding rules appearing the columns and the four different non-response treatments in the rows.^[
Complete details on these simulations are provided in the online Supplementary Materials.]

The upper left pane simulates the combination of the mixed coding rule and the treatment of survey non-responses as indicating a lack of support for democracy matches, the same combination employed in @Claassen2020b.
The initial drop and slow recovery in the mean of these simulations of public democratic support was the evidence presented for the "thermostatic response of public opinion" and the claim that there is "little evidence that democracy generates its own demand" [@Claassen2020b, 48].
But the other panes, with a few exceptions, show smaller dips, quicker recoveries, and continued increases; these findings lead to very different conclusions.
Indeed, most of the analyses---eleven of sixteen---show statistically significant increases in democratic support within three decades, the sort of generational change predicted by the classic theory since @Lipset1959.
None of them show statistically significant declines that would lend credence to the argument that democratic support responds thermostatically.


## Discussion {.unlisted .unnumbered}

@Claassen2020b claims that a thermostatic response to democratic change is required to explain "democracy’s fading allure" among publics.
It is far from clear that support for democracy is indeed declining [for recent evidence against, see, e.g., @Wuttke2022].
But if it is, the fragility of the results presented in @Claassen2020b to alternate reasonable coding decisions and non-response treatments suggests that other explanations are worth revisiting.
One possibility not considered in the specification adopted in @Claassen2020b is that rising economic inequality undermines support for democratic institutions [@Magalhaes2014, 84; see also, e.g., @Solt2012].
This and other theories of democratic support warrant further investigation.

Regardless, the fragility of the findings in @Claassen2020b that is documented above has profound implications.
The thermostatic theory of democratic support presented in @Claassen2020b suggests that would-be autocrats working to erode democracy constitute little cause for concern: "should elected leaders start dismantling democratic institutions and rights, public mood is likely to swing rapidly toward democracy again, providing something of an obstacle to democratic backsliding" [@Claassen2020b, 51].
To the contrary, that theory implies, it is those who attempt to deepen democracy who pose the threat, as "extending democratic rights and legal protections to minorities" will trigger a public backlash against democracy [@Claassen2020b, 51].
If the theory is true, the appropriate positions for those who value democracy are, in short, complacency and satisfaction with the status quo.

That the evidence for such a perverse outcome depends entirely on researcher degrees of freedom indicates that defending democracy requires taking a very different course.
Alternate combinations of plausible coding rules and non-response treatments yield no support for the thermostatic theory.
In fact, most show that increases in democracy work, over time, to increase the public's support, just as suggested by the classic theory that has been advanced in scholarship over more than a half century.
These findings imply that those who want to advance democracy should do just that: they should push ahead.

## Acknowledgments {.unlisted .unnumbered}


Yue Hu acknowledges the support of the National Natural Science Foundation of China [72374116].
Frederick Solt acknowledges the support of the Office of the Provost, University of Iowa.
The three authors thank the editors and anonymous reviewers for their helpful comments.

## References {.unlisted .unnumbered}

::: {#refs}
:::

\pagebreak

\newpage
\appendix

# Online Supplementary Materials {.unlisted .unnumbered}
\captionsetup[figure]{list=yes}
\setcounter{page}{1}
\renewcommand{\thepage}{SI-\arabic{page}}
\setcounter{figure}{0}
\renewcommand{\thefigure}{SI.\arabic{figure}}
\setcounter{table}{0}
\renewcommand{\thetable}{SI.\arabic{table}}
\tableofcontents
<!-- \listoftables -->
<!-- \listoffigures -->

\newpage
# Tabular Results

```{r num-labelAPSR}
names_coef <- c(
  "(Intercept)" = "(Intercept)",
  "plm::lag(SupDem_trim, 1:2)1" = "Democratic Mood (t-1)",
  "plm::lag(SupDem_trim, 1:2)2" = "Democratic Mood (t-2)",
  "diff(Libdem_z, lag = 1)" = "Liberal Democracy (Difference)",
  "Libdem_z" = "Liberal Democracy (Difference)",
  "plm::lag(Libdem_z, 1)" = "Liberal Democracy (t-1)",
  "diff(Polyarchy_z, lag = 1)" = "Electoral Democracy (Difference)",
  "Polyarchy_z" = "Electoral Democracy (Difference)",
  "plm::lag(Polyarchy_z, 1)" = "Electoral Democracy (t-1)",
  "diff(Liberal_z, lag = 1)" = "Minoritarian Democracy (Difference)",
  "Liberal_z" = "Minoritarian Democracy (Difference)",
  "plm::lag(Liberal_z, 1)" = "Minoritarian Democracy (t-1)",
  "diff(lnGDP_imp, lag = 1)" = "Log GDP Per Capita (Difference)",
  "lnGDP_imp" = "Log GDP Per Capita (Difference)",
  "plm::lag(lnGDP_imp, 1)" = "Log GDP (t-1)",
  "diff(Corrup_TI_z, lag = 1)" = "Corruption (Difference)",
  "Corrup_TI_z" = "Corruption (Difference)",
  "plm::lag(Corrup_TI_z, 1)" = "Corruption (t-1)"
)

names_gof <- tibble::tribble(
  ~raw,           ~clean, ~fmt,
  "r.squared", "$R^2$",    3,
  "rmse",  "RMSE",    2
)
```

```{r m11tab1, results='asis'}
m1_1_modelsummary <- m1_1_list %>%
  unlist(recursive = FALSE)

modelsummary(m1_1_list[1:2],
             shape = "rbind",
             stars = TRUE,
             coef_map = names_coef,
             gof_map = names_gof, 
             coef_omit = "Int",
             output = "latex", 
             title = "Predicting the Change in Public Democratic Support") %>% 
  kableExtra::add_header_above(c(" " = 1, "Non-Response Treatment" = 4)) %>% 
  kableExtra::kable_styling(font_size = 7) %>% 
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

```{r m111tab2, results='asis'}
modelsummary(m1_1_list[3:4],
             shape = "rbind",
             stars = TRUE,
             coef_map = names_coef,
             gof_map = names_gof, 
             coef_omit = "Int",
             output = "latex") %>% 
  kableExtra::add_header_above(c(" " = 1, "Non-Response Treatment" = 4)) %>% 
  kableExtra::kable_styling(font_size = 7) %>% 
  kableExtra::kable_styling(latex_options = "HOLD_position")
```
<!-- \elandscape -->

\newpage
# Calculating Effects Via Simulation

To simulate the effects of changes in democracy on public support for democracy in the error-correction models [see @Williams2012], we follow the same procedure employed in @Claassen2020b [, 48-50].
That is, for each of our sixteen sets of results, we first set all of the independent variables to the same moderate values used in @Claassen2020b and ran the model for 200 years, which is a sufficient time for the system of equations to stabilize.
Next, the value of democracy was increased one standard deviation, from its previous value of a half standard deviation below its mean to a half standard deviation above its mean.
Finally, the system of equations was allowed to run for thirty additional years.
These thirty years of simulated results are shown in the sixteen panes of Figure \@ref(fig:simulationPlot) in the text.
Again following @Claassen2020b [, Supplementary Information 3] and @Claassen2020c, the uncertainty in the model was captured by taking 10,000 draws from a multivariate normal distribution with the expectation that the vector of model coefficients and variance constitute the robust covariance matrix, $\tilde{\Theta} ∼ MVN(\Theta, \Sigma)$, and adding the noise estimated in the regression standard error, $\tilde{Y}_i ∼ N(X_k \tilde{\Theta}_{ki}, \sigma)$.
To get first differences, the mean value of $\tilde{Y}_i$ in the year before the increase in democracy ($t = -1$) was subtracted from each $\tilde{Y}_i$, and the 0.025 and 0.975 quantiles of the first difference were used as its lower and upper confidence bounds.
